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ABSTRACT 

TeV-blazars potentially heat the intergalactic medium (IGM) as their gamma rays interact with 
photons of the extragalactic background light to produce electron-positron pairs, which lose their 
kinetic energy to the surrounding medium through plasma instabilities. This results in a heating 
mechanism that is only weakly sensitive to the local density, and therefore approximately spatially 
uniform, naturally producing an inverted temperature-density relation in underdense regions. In this 
paper we go beyond the approximation of uniform heating and quantify the heating rate fluctuations 
due to the clustered distribution of blazars and how this impacts on the thermal history of the IGM. 

We analytically compute a filtering function that relates the heating rate fluctuations to the underlying 
dark matter density field. We implement it in the cosmological code GADGET-3 and perform large 
scale simulations to determine the impact of inhomogeneous heating. We show that, because of blazar 
clustering, blazar heating is inhomogeneous for z > 2. At high redshift, the temperature-density 
relation shows an important scatter and presents a low temperature envelope of unheated regions, 
in particular at low densities and within voids. However, the median temperature of the IGM is 
close to that in the uniform case, albeit slightly lower at low redshift. We find that blazar heating is 
more complex than initially assumed and that the temperature-density relation is not unique. Our 
analytic model for the heating rate fluctuations couples well with large scale simulations and provides 
a cost-effective alternative to subgrid models. 

Subject headings: BL Lacertae objects: general, cosmology: theory, gamma-rays: general, intergalactic 
medium, large-scale structure of universe 


1. INTRODUCTION 
The intergalactic medium constit utes t he b ulk of the 


baryons forming the cosmic web (Bond et al. 

1996) as 

well as the reservoir of baryons available for tr 

le forma- 


tion of galaxies and clusters (Rauch et al. 1997). Ob¬ 


servations of metal-enriched gas link it s evolution to the 


formation of galaxies and st ars (see e.g. Simionescu et al 


(2009); Werner et al. (2010)). Therefore, the thermal his 


tory of the 1G1V1 plays a central role in determining the 
development of structure in the visible universe. 

The temperature of the IGM is mostly set by photoion¬ 
ization of hydrogen and helium, competing with adia¬ 
batic cooling. As a result, the universe is slowly cooling 
once reionization is completed. The IGM temperature 
is expected to increase with density because denser re¬ 
gions experience a lower amount of adiabatic cooling, as 
well as more recombination-induced photoheating. Gas 
in the IGM that has not been shock heated provides a 
lower envelope to the temperature-density distribution. 
Observations of Lyman a absorption lines with 4 ^ z ^ 2 
show that lower column density gas exhibits the lowest 
line widths, which is commonly interpreted as an indi- 


gas (Kirkman & Tytler 1997! 

Schaye et al. 2000 Ricotti 

et al. 2000; Rudie et al. 20E 

ij). A temporary flattening 
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of the slope of the inferred temperature-density relation 
around z ~ 3 indicates an additional source of heating 
at that redshift. 

The latter is expected due to Hell reionization (e.g 


Tittl ey fc Meiksmj 20071 IMcQuinn et al.| 2009| |Com 
postella et al.||2013[ Puchwein et al. 2014). Observations 
suggest that the tail end of He II reiomzation occurred 


between z ~ 3.2 and 2.7 (Kriss et al. 2001 Syphers & 


Shull 2014|). However, the bulk of He 11 has likely ion¬ 
ized e arlier, potentially as early as z > 4 (jWorseck et al. 


2014). Its patchiness is probably related to the rarity 
of Ionizing sources and the inhomogeneity of the IGM. 
Including full radiative transfer, simulations show en¬ 
hanced scatte r in the temperature- density distribution 
around z = 3 ( McQuinn et al.|[2009 ). Whether reioniza¬ 
tion ma y result in an inverted te mperature-density distri- 
buti on ([M eiksin & Tittley [2012 ) is not firmly established 
yet (Compostella et al. pT3l ). 

On top of that, recent measurements found that un- 


derdense regions may be warme r than predicted (Viel 
Bolton et al. 2008), as well as to exhibit 


et al. 2009 


unexp 

2014]), which is harder to explain solely by ionization "of 
He 11. Although an inverted temperature-density relation 
in un derdense regions has not been firmly established yet 
(Bolton et al.] |2014 ), these observations suggest that the 
thermal history of the IGM may be more complex than 


ectedly high temperatures at z < 2 (Boera et al. 
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i nitially assu med. 


Broderick et al. (2012) recently suggested a comple¬ 


mentary heating mechanism, through TeV blazars. TeV 
blazars are active galactic nuclei (AGN) emitting very 
high energy gamma rays (E > 100 GeV). They be¬ 
long to the radio-loud subgroup of AGN, with the rel¬ 
ativistic jet being pointed towards us. About 50 of 
these sources have been significantly discovered so far 
(http://tevcat.uchicago.edu/) by ground based Cerenkov 
telescopes such as MAGIC, H.E.S.S. and VERITAS. 
Those pointed observations have only skimmed the sur¬ 
face of a much larger population that manifest them¬ 
selves as hard-spectra gamma-ray blazars as observ ed 
with the sp ace-based Fermi -LAT telescope (Broderick 
et al.||2014a|). 

The universe is mainly opaque to very high energy 
gamma rays; they interact with the extragalactic back¬ 
ground light (EBL) producing electron/positro n pairs 


(Gould & Schreder 1967 Stecker et al. 1992). It is 


commonly assumed that the electron/positron pairs in¬ 
verse Compton scatter photons of the cosmic microwave 
background, resulting in a distribution of photons with 
energies between 0.1 and 100 GeV. Such an emission 
compon ent towar ds TeV blazars has not been observed 
so far (Aleksic et al. 2010|), despite ex tended searches 
(H. E. S. S. Collaboration et al. 12014). One solution 
would be pair deflection due to the intergalactic mag¬ 
netic field, thus lowering the surfac e brightness of the 
formed pair halo at GeV energies dDurrer & Neronov 


2013 Vovk et al. 


Dermer et al. 2011 j" 


However, the cascaded GeV emission would still con¬ 
tribute to the extragalactic gamma-ray background 
(EGRB). Fermi -LAT was able to resolve more sources 
than its predecessor EGRET, thereby limiting the 
isotropic component of the unresolved EGRB and 
severely constraining the redshift evolution of hard 


Venters 2010 Murase et al. 

2012 

Inoue & Ioka 

2012). 

As a result of these, it is now weJ 
this picture, the co-moving numbe 
ray blazars, and by extension the ^ 
ther be constant or decreasing witf 

il established that in 
;r density of gamma- 
re V blazars, must ei- 

i redshift dKneiske & 

Mannheim 2008; Venters 2010; Abazajian et al. 2011; jin- 


oue v Ioka 

classes of AGNs* specifically, and all other tracers of 
the cosmological history of accretion onto galactic ha¬ 
los generally (e.g., star formation). Furthermore, af¬ 
ter careful modeling of gamma-ray and X-ray survey 
selection effects there appears to be no evidence for a 
significantly different evolution of blazars in compari- 
son to their radio-loud an alogues such as radio galaxies 


(Giommi et al. 2012, 2013). Together this casts doubt 
on the inverse Compton cascade picture and provides 
circumstantial evidence that pair beams instead transfer 
their en ergy directly to the IGM through plasma insta- 


regim e extends to the non-linear regime, | Chang et al. 
(2012) show it is responsible for increasing the temper- 
ature of the IGM by almost a factor 10 in low den¬ 
sity regions. Wh il e this assumption is still d e bated (see 


Miniati fc Elyiv (|2013); [Sironi fc Giannios 


. (2 on 

also Schlickeiser et al. ~j 20l!H ~ 2012p ; Chang et al. ~ p014j )) 


but 


throughout all this paper we assume plasma mstabili- 
ties are the dominant mechanism for cooling of the pair 
beams. 

Including TeV blazar heating in the thermal history 
of the IGM, Chang et al. (2012) were able to reproduce 
the inver ted temperature- density relation for low density 


regions. Pfrommer et al. (2012) found that TeV blazar 
heating is capable of creating a redshift dependent en¬ 
tropy floor in clusters and galaxies, thus suppressing the 
formation of dwarf galaxies after the peak of blazar activ¬ 
ity at redshift z ~ 2 and potentially providing an expla¬ 
nation to the “missing sa tellite problem ” and the “miss¬ 
ing void dwarf problem” (|Kravtsov[2010|. Implementing 
uniform blazar heating, ne! with a redshift-dependent 
but spatially homogeneous energy deposition rate per 
unit volume, in a cosmol ogical hydro dynamical simula¬ 
tion of galaxy formation, Puchwein et al. (2012) find ex¬ 
cellent agreement with the one and two-pomt statistics 
of the Lyman a forest, which is the main observational 
tracer of low density regions in the universe. 

However, the low thermal broadening of certain Ly- 
man a lines indicat es the presence of cold gas at z = 2.4 
(Rudie et al. 2012), which suggests TeV blazar heating 
does not uniformly heat the whole universe. It is natural 
to expect a larger TeV flux close to higher density re¬ 
gions where visible structures form. Conversely, in large 
under dense regions, far from massive black holes, heat¬ 
ing is probably much lower. The goal of this paper is 
to go beyond the hypothesis of uniform heating and to 
link TeV blazar heating to the underlying clustered den¬ 
sity field and take into account the bias of sources. This 
will lead to a more heterogeneous heating pattern and 
account for unheated regions while keeping the overall 
impact of blazar heating. 

Self-consistently studying the evolution of the IGM 
from first principles involves modeling both the forma¬ 
tion and evolution of galaxies at the largest scales of the 
universe. As this is still far beyond reach of current com¬ 
puters, we have determined a filter function which relates 
the heating fluctu ations to the dark mat ter (DM]_struc 
ture sim i larly to [Pritchard fc Fu rlanettol (|2007|) ; [Barkana 


& Loeb (2005); Pontzen 12014) 


Based on the hierarchi¬ 
cal structureT’ormation inliACDM universe, it naturally 
selects the relevant length scales for TeV blazar heating 
(§2). We have implemented it in large scale cosmological 
simulations (§3) in order to focus on the equation of state 
and thermal evolution of the IGM (§4). We then discuss 
how inhomogeneous heating could reconcile different ob- 


bilities (Broderick et al. 

2012 

Schlickeiser et al. 

2012 


2013; 

Chang et al. 2014 

). This naturally explains the 2 . DETERMINING THE WINDOW FUNCTION 


evolution in agreement with that of quasars (|Broderick 
et al. 1|2014a|b[ ). 


2 . 1 . 


The pairs constitute a dilute, ultrarelativistic beam, 
which is subject to several pl asma instab i lities , from 
which the “oblique” instability (Bret et al. 2004) is the 
most powerful. Assuming its efficiency in the linear 


One zone models (Chang et al. 
2012) and numerical simulations - 


2012 j iPfrommer et al . 

( jPuchwein et al. 
of blazar heating on the IGM assume that the heating 
is uniform. Because the heating rate depends on the lo¬ 
cal density of EBL and TeV photons, the assumption of 
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uniform heating implies that the distributions of EBL 
and TeV photons are uniform. The EBL photon density 


in voids is only 2% less than the average value (Furniss 
et al. 2015). However, for TeV photons the mean free 


path compares with the separation between TeV sources 
for z ^ 1, so the spatial fluctuations in the heating rate 
are likely nontrivial. Moreover, the sources of TeV pho¬ 
tons tend to be clustered and so the IGM near these 
clustered regions will get an increased flux of photons 
in comparison to low-density regions because of the in¬ 
creased number of sources and the 1/r 2 flux dilution with 
increasing distance r from the source. 

Our goal is to include a more realistic model for heat¬ 
ing due to TeV blazars in numerical simulations. To 
properly calculate the heating fluctuations due to TeV 
blazar heating, the formation and evolution of accreting 
supermassive black holes must be modeled in a full self- 
consistent cosmological simulation. In addition, the TeV 
radiation from these systems must be ray-traced through 
the simulation volume. Such a task is computationally 
intractable. As a result, we have elected to model this 
TeV blazar heating in a statistical manner. 

We assume that TeV blazars are associated with galax¬ 
ies or alternatively quasars and that they roughly emit 
over 47r steradian. The latter assumption remains valid 
as long as the duty cycle of blazars is small enough so that 
jets point in all directions over cosmological timescales 
and the number of TeV blazars is large enough such that 
every spot in the universe is illuminated by at l east a few 


TeV blazars, which is the case for z < 6 (Chang et al 
2012). The heating rate at a given point x is determined 


2.2. Window function for TeV blazar heating 

To determine the heating rate fluctuations we express 
the TeV emissivity in Eq. (fTl) as a mean value and a first 
order correction. The heating rate at a given point is 
set by the received TeV flux from all the sources within 
a certain radius. We assume that the pairs lose their 
energy to the IGM at the point wher e the y are created. 
As sta ted in Broderick et al. (2012) and Chang et al. 
(2014), this is a reasonable assumption as the plasma 
instability length scales are many orders of magnitude 
smaller than the mean free path of these TeV photons. 

The TeV gamma rays are emitted by accreting super- 
massive black holes at centers of galaxies, which cluster 
in overdense regions. Matter is tightly coupled to the un¬ 
derlying dark matter, the evolution of which is straight¬ 
forward to model analytically within the linear approxi¬ 
mation. The linear approximation is valid as long as the 
overdensity is small, which is true in the early universe 
and then breaks down at small scales as very dense struc¬ 
tures form. Our computation takes into accountjthe bias 
betwe en baryonic matter and dark matter (Mo & White 
1996), as we detail below. 


by the received TeV flux from all the sources 


To model cosmic distances, Eq. 0 is integrated in red- 
shift space and we take into account the resulting energy 
loss for the TeV photons as well as first order correc¬ 
tions due to a clustered blazar distribution and proper 
motions of the sources within the Hubble flow (Kaiser 
1987). We integrate over the energy distribution of the 
TeV-emission. Following the detailed derivation in the 
Appendix, we obtain the window function H/#(&, z) as a 
function of comoving wavelength 




,£(r + x, z ) _ T 
dr n 77\ e ’ W 

0 ^PP \ z ) 


where Q and £ are the heating rate density and emissivity 
of the sources (both in units of energy per unit time and 
per unit volume), D vv> is the mean free path for pair 
production on the EBL, r is the associated optical depth 
along the line of sight and z' is the redshift corresponding 
to the distance |r'|. 

One can then express the resulting heating rate as a 
mean value Q and a first order correction 5h- 


Q(x,z) = Q(x,z)[1 + S h (x,z)\. (2) 

The method is based on a Taylor expansion of the quanti¬ 
ties describing the TeV sources and keeping only the first 
order corrections. Transforming to Fourier space yields 
the fluctuation amplitude 

S H (k,z) = W H (k,z)5(k,z), (3) 

where Wh is defined as the window function and maps 
the Fourier transform of the dark matter overdensity, S, 
to the Fourier transform of the heating fluctuations, 5h- 
This naturally yields the relevant length scale for heating 
rate fluctuations. The detailed and pedagogical exposi¬ 
tion of this method (in the Newtonian limit, in an ex¬ 
panding universe, and additionally accounting for a clus¬ 
tered source distribution) is given in the Appendices [A] 
to [C] In the following section, we present the general 
method and highlight the underlying hypotheses of our 
work. 


W H (k,z) = 


i r E 

N.Ie„ 

DV) 

D(z) 


dE 


D PP (E , z) 


i; 


dz' 


K z> ) + 0 jo(kr) ~ yj 2 (hr) 


£ E '{z')e- T 
(1 + z')H(z') 

( 4 ) 


with 


N = 


/ 

Je, 


dE 


D pp {E,z) 


i; 


dz' 


£ E '{z’)e T 
(1 + z')H{z'Y 


( 5 ) 


Here, f is the comoving distance between the source and 
heated regio n, D denotes the linear growth factor defined 
in Eq. (B8), / = dlog£/dloga, E is the energy of the 
received TeV photon, E' its initial energy, and £e the 
comoving blazar spectral luminosity density, b is the Eu- 
lerian bias and jo and are spherical Bessel functions. 

The blazar spectral luminosity density (in units of en¬ 
ergy, per unit time, per volume, per energy) is deter¬ 
mined by both the intrinsic source spectra and the lu¬ 
minosity function o f TeV blazars. Here we adopt the 
model described in Broderick et al. (2014a), in which 
the spectra are assumed to be well described by a set 
of broken power laws and a luminosity function that de¬ 
pends on redshift, luminosity from 100 GeV to 10 TeV 
(Lxev)j and the low energy spectral index (T/). The 
assumed source spectral luminosity density is given by 
Te (AxeV? T/) = LE(Vi)L Te v where 


f (r \ = _M_ 

M 11 (E/E b ) r ‘ + (E/E b ) r » ’ 


( 6 ) 






























4 


where the normalization /o is 


fo.inZ^BdE = 1, E b = 1 TeV, 


chosen such that 

.GeV ^ , and = 3. The 

luminosity function is based on the spectral properties 
of the Fermi population of hard gamma-ray blazars a nd 
the qu asar luminosity function reported by |Hopkins et al. 

2007). In terms of these, the blazar spectral luminosity 


ensity is 


/ 


£e(z) = I d logio ^TeV dTI L E (L Te v,ri)<j)(z, I/TeV, T/) , 

(7) 

where is the comoving analog of Eq. 6 in 

Broderick et al. (2014a). This may be simplified further 


using the separability of 0 in T/, i.e., 

T TeV , Ti) = cj)(z, L TeV ) X (Yi ), ( 8 ) 

where f dTix(Fi) = 1 . As a result, 


where 


&e(z) = M z ) > 

l e ) = J dTt l e ( r,)x(rO 


which is now a function of energy alone, and 
k(z)= / d\og w L TeV 4>(z,L TeV ) 


(9) 


( 10 ) 


( 11 ) 


is the TeV luminosity density of TeV blazars. 

In practice, for the T eV blazar luminosity function in 


Broderick et al. (2014a) \ L>e] is well fit over the energy 

range of interest, 100 GeV to 10 TeV, by Eq. J6|) with 
Ti = 1.80 and E & = 1.057, with an L\ norm of0.0048. 
This spectral shape peaks near 1 TeV and thus sets a 
characteristic 7 -ray energy. 

Sinc e the T eV blazar luminosity function in ! 


et al. (2014a]). _is simply the quasar luminosity function in 
Hopkins et al. (2007) rescaled in energy and luminosity, 
the comoving luminosity densities are also proportional. 
That is, 


k(z) = CA Q (z), 


( 12 ) 


with £ = 2.1 x 10 3 and the comoving quasar luminosity 
density A q. The co rresponding co movi ng TeV luminos¬ 
ity density, found by Chang et al. ( 2012 ) b ased on a fit to 
the co moving quasar luminosity density in|Hopkins et al. 
p007| , is given by 


where Aq = (1.7 — 4.8) x 10 36 erg s 1 cm 


A(z) = Ao10 1,18 ^ _0,081222_0,182 ^ 3+0,051524_0,0041825 

(13) 

i -3 the blazar 

luminosity density at the current epoch. 

The optical depth is set by the physical mean free path 
D pp (E',z) of TeV photons before they interact with an 
EBL photon to produce an electron-positron pair. Its 
redshift evolution is set by the density of the EBL. De¬ 
spite careful studies that aim at constraining it, there 
remain uncertainties in the star formation history of the 
universe as well as its metallicity and dust contents (see, 


e.g. 

ing 


Franceschini et al.|2 008 Stecker et al.|2006). Follow- 
Uhang et al. ( | 2012 | ) we use a prescription 


D pp (E',z') = 35 


E' 


1 TeV 


1 + z 


-e 


Mpc, (14) 


where £ = 4.5 for z < 1 and £ = 0 for z > 1 (Kneiske 
et al. 2004 Neronov & Semikoz] 2009) and the mean tree 
path is expressed in physical units. The proper mean 
free path is constant for 2 > 1 , however, the comoving 
mean free path increases for increasing redshift. 

The fluctuations in TeV flux are related to the distri¬ 
bution of blazars, which is biased with respect to the 
distribution of dark matter halos. Luminous structures 
such as galaxies preferentially populate the high peaks of 
the dark matter density distribution. The square of the 
bias b of a given structure is the ratio between its power 
spectrum to the power spectrum of the dark matter halos 
and it is stronger for more massi ve objects, such as the 
host galaxies of quasars (see, e.g. Cooray & Sheth 2002 , 
for a review). 

Due to the small number of sources and TeV photon 
absorption, there is no observation of TeV blazar bias. 
Therefore, we use a model for quasar bias as well as a 
model for galaxy bias to illustrate the impact of the un¬ 
certainty on the value of the bias. Since AGN are gener¬ 
ally much more biased than galaxies the latter provides 
a robust lower bound. However, TeV blazars may have a 
stronger bias than quasars, as radio-loud AGN are gener¬ 
ally found in a more clustered environment than q uasars 


(Mandelbaum et al. 2009 Simpson et al. 2012 ). Our 


m odel for quasar bias f o r z > 1 is ba s ed on a fit of data 


by_ Croom et al. 
( 2007| ). We use 


(2005); Myers et al. (2007); Shen et al 


r(z) = 10 1 


,0.272-0.04 


(15) 


which is very similar to Papageorgiou et al. ( 2012 )_. Our 
model for galaxy bias is based on a ti t of data by Mari- 
noni e t al. (2005); Steidel et al. (1998); Kashikawa et al. 
We use 


^galaxy (-^) — 10 


,0.1742 


(16) 


Figure [l] shows the corresponding evolution of quasar and 
galax y bias as a function of redshift. Using the same 
data, Basilakos et al. ([2008) find that the quasar and 
galaxy bias can be fitted using a dark matter halo model 
of 10 13 h~ 1 M Q and 10 ^ 12 h~ 1 AE 1 _ respectively. 

Substituting Eqs. Q and (14) into Eq. Q then gives 


the complete window function for TeV blazar heating. 
Figure [2] shows the filter for z = 1,2 and 4 for a model 
with galaxy and quasar bias. We have computed the win¬ 
dow function using an embedded Runge-Kutta method 
which is able to capture the fast variation of the Bessel 
functions at large wave numbers while decreasing com¬ 
puting time at smaller wave numbers. The window func¬ 
tion is integrated up to z = 8 , which ensures that all the 
sources are taken into account. Changes in maximum 
redshift of the integration have no discernible impact on 
the window function. 

The window function describes how density fluctua¬ 
tions translate into heating fluctuations. The k _1 slope 
for large wavenumbers results from the r ~ 2 decrease in 
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Fig. 1.— Galaxy (red line) and quasar (blue line) bias models 
used in our simulations. 


10 3 
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Fig. 2.— Window function for TeV blazar heating from z = 1 
(solid lines), z = 2 (dashed lines) and z = 4 (dotted lines) for the 
galaxy bias model (red) and the quasar bias model (blue). The 
vertical dashed lines indicate the minimal and maximal comoving 
wavenumber modeled in the simulation. 


the TeV-flux. The position of the break is set by a com¬ 
bination between the mean free path and the redshift 
evolution of the blazar luminosity density. One might 
expect the comoving mean free path, which increases for 
increasing redshift at z > 1, would move the break to 
smaller wavenumbers at higher z. However, by integrat¬ 
ing over the energy of the photons and including lower 
energy photons that have longer mean free paths, the 
effect of the change in the mean free path is mitigated. 
This allows the rapidly decreasing blazar luminosity den¬ 
sity for z > 2 to move the break to larger wavenumber: 
the luminosity density drops too steeply so that more 
distant sources within the photon mean free path do not 
contribute as much as sources more close by. This lat¬ 
ter effect moves the break to larger wavenumber for in¬ 
creasing redshift when lower energy photons with longer 
mean free paths are accounted for. This implicitly as¬ 


sumes that these low energy photons would contribute 
to the local heating rate via pair production and plasma 
instabilities. A more rigorous calculation that accounts 
for the efficiency of blazar heating would vary with pho- 
ton_energy and redshift (see for instance Chang et al. 
2014 for the effect of nonlinear Landau damping). How- 
ever, this efficiency is not yet fully understood and so 
we assume in this paper that it as a function of pho¬ 
ton energy is unity. As the normalization is set by the 
bias, the quasar bias model has more power at all scales 
and density fluctuations will lead to more enhanced heat¬ 
ing. In both models, at high redshift, most of the power 
resides on large scales, where blazar heating traces the 
density fluctuations. At smaller scales, density fluctua¬ 
tions have no impact and blazar heating is uniform. At 
the current epoch, TeV blazar heating is close to uniform 
because there is power only at the largest scales (above 
100 M pc) wh ere t he universe is essentially uniform (see 
Clowes et ah] (2013) and references therein). 

In both models, the window function remains posi¬ 
tive at all scales; thus, underdense regions are the only 
areas where a lower than average heating rate is ex¬ 
pected. Heating at rates larger than the average is pos¬ 
sible for large-scale overdensities but also on small scales 
for highly non-linear objects with S 1. We caution 
that our simplified approach with a linear bias descrip¬ 
tion starts to break down there and would have to be 
augmented with a non-linear description of bias. How¬ 
ever, the gas in these regions experiences shock heat¬ 
ing in addition to photoheating such that the influence 
of blazar heating becomes negligible there. Moreover, 
those densities have little influence on the statistics of 
the high-redshift Lyman a forest and thus shall not be 
the subject of the present work. After these first analytic 
estimates, we include the window function to model TeV 
blazar heating in cosmological simulations. 


3. NUMERICAL METHOD 
3.1. Cosmological simulations 

We perform simulations with the smoothed particle hy¬ 
drodynamics (SPH) code GADGET-3, an upgraded ver 


sion of the publicly available GADGET-2 code (Springel 
2005). The code solves the gravitational evolution of 


both dark matter and gas particles following a TreePM 
N-body method. The hydro dynamical evolution of the 
gas is m odeled u sing an entropy conserving scheme 
(Springel & Hernquistj2002). 

'Th e cos mological m odel is based on the WMAP 7-year 
data (Komatsu et al. 2011): VLm = 0.272, = 0.728, 

Qb = 0.0465, h = 0 704 and erg = 0.809. The initial 
conditions were evolved from z = 100 until z = 1 in 
boxes with comoving side length of 100 h -1 Mpc and 
periodic boundary conditions. We use N = 2 x 512 3 par¬ 
ticles, which gives a mass of m gas = 3.8 x lO 6 h _1 M 0 and 
ttidm = 1.8 x 10 7 h~ 1 M® for baryonic and dark matter 
particles, respectively. We used a comoving gravitational 
softening length of 7.8 h~ x kpc. We checked that the res¬ 
olution has a barely discernible impact on the results of 
our simulations by performing test simulations with a 
comoving side length of 50 h~ x Mpc at different resolu¬ 
tions, up to 2 x 512 3 . The spread of the temperature 
in the low density regions increases with increasing res¬ 
olution. However, this effect is approaching saturation 
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at a resolution of N = 2 x 256 3 . We find less than 4% 
difference in the median temperature and less than 10% 
difference in its root mean square in low density regions 
between N = 2 x 256 3 and 2 x 512 3 simulations, indicat¬ 
ing that our chosen resolution is sufficient to accurately 
capture the temperature-density distribution. 

As we are only interested in the low density intergalac- 
tic medium, we use a simplified model for star formation 
which significantly speeds up the simulations. In this 
model, gas particles with £ g as > 1 000 and T < 10 5 are 


directly converted into stars (Viel et al.|2004). Although 
it results in unrealistic galaxy properties, this approxima¬ 
tion does not affect regions with S gas <10. Local black 
hole feedback is not included. Photoheating is set by ion¬ 
ization equilibrium of H, He I and He II in the presence of 
an external UV fi eld, whic h is p arameterized according 
to 


_para: 

Faucher-Giguere et al.| (2009). As our version of the 


GADGET-3 code assumes ionization equilibrium when 
computing photoheating rates, the heating is rather in¬ 
efficient during reioniz ation where this assum ption is not 
well satisfied (se e e.g. | Puchwein et al. 2014). Following 
Puchwein et al.| (|2012|) we thus include the equivalent 
heat input by hand at redshift z = 10. Our simulation 
does not include radiative transfer effects on the pho- 
toionization of He II. 

The size of the box is set to model the heating pertur¬ 
bations on the scales determined by the window function 
in Figure |2j We want to model a representative cosmic 
sample and probe distances beyond the mean free path 
of the TeV photons, which is of order of 70 (comoving) 
Mpc at z = 1 (but 105 Mpc at z m 2). To confirm that 
100 h~ x Mpc is a satisfactory size to model all the signif¬ 
icant length scale, we performed a Tbox = 200 h~ x Mpc 
simulation with 512 3 particles. We compared the result¬ 
ing temperature distribution function with the 100 h~ x 
Mpc box, with the same mass resolution. The mean tem¬ 
perature varies by less than 2% for regions with S gas ^ 0 
at all redshifts. The standard deviation varies by 10% 
at z = 3 and is below 5% at z = 1. This is consistent 
with studies showing that the one-and two-point statis¬ 
tics of the Lyman a forest are well captured with ~ 50 
h -1 boxes ( Regan et al.||2007 Bolton fc Becker||2009 ). 

3.2. Including the TeV blazar heating fluctuations 

We model the impact of the fluctuations in TeV blazar 
heating on the thermodynamics of the IGM. For every 
gas particle, the blazar heating is set by the mean value 
plus some correct ion depending on the lo cal density field 


(Eq. ([2|)). As in Puchwein et al. (|2012|), we adopt the 
mean "heating rate computed by [Chang et al. (2012) 
(Eq. ([ 9 ])) . We focus on the model with “intermediate 7 ^ 
values tor the blazar heating. 

Modeling fluctuations by implementing a filtering func¬ 
tion in a large scale simulation is a new method. The 
computation of the fluctuations is done in Fourier space 
and is inspired by the numerical particle-mesh (PM) al¬ 
gorithm used to solve the long-range gravitational force. 
The Fourier transforms are performed with the parallel 
extension of the Fast Fourier Transform Library. The 
first step is map the particles onto a mesh, which is 
done with a clouds-in-cells algorithm (Hockney & East- 


wood 1981). We determine the Fourier transform of the 
DM density field. Then the density field is multiplied 
by the window function performing a bilinear interpo¬ 


lation of tabulated values for certain values of redshift 
and wavenumber. In this paper we use 21 equally spaced 
redshift bins from z = 5, where blazar heating turns 
on in our model, until the end of the simulation. We 
use 128 logarithmically equally spaced wavenumber bins. 
Similarly to the long-range gravitational PM force, we 
deconvolve for the clouds-in-cells kernel by dividing by 
sine 2 (k x L/ 27V)sinc 2 (k y L/2 7V)sinc 2 ( k z L /2 N) . We then 
perform the inverse Fourier transform and renormalize. 
The last step is necessary to remap the results onto the 
gas particles. As our filtering technique may result, in 
rare cases, in unphysical cooling (Q < 0 ), we set a lower 
limit of S H = — 1 . 

4. RESULTS 

Figure [3] shows the heating rate fluctuations in the mid¬ 
plane of the L x = 100h -1 Mpc simulation for z = 3 and 
1 for both the galaxy and quasar bias model. Both simu¬ 
lation models are started from identical initial conditions 
such that any visible difference is solely due to the differ¬ 
ent bias assumptions in the blazar heating model. The 
corresponding density field is shown in the upper column 
and shows increasing structure formation as the redshift 
decreases. The heating map has a linear scale while the 
density scale is logarithmic. The heating rate fluctua¬ 
tions are, on average, much smaller and more spatially 
homogeneous than the density fluctuations. This is be¬ 
cause the window function filters out small scales, which 
correspond to collapsed regions, where density fluctua¬ 
tions are the highest. To the zeroth order, one can thus 
consider TeV blazar heating to be uniform, as was as¬ 
sumed in Chang et al. ( 2012 ). 

Additional heating (i.e. 6h > 0 ) occurs around clus¬ 
tered regions. This is expected, as the window function 
translates large scale density fluctuations into heating 
rate fluctuations. Conversely, under dense regions, such 
as the one around {x = 60, y = 70} (for z = 3) display 
heating below average as they are isolated from sources, 
and their flux decreases as e~ T r~ 2 . This is most obvious 
at high redshift, and for the quasar model, where bias 
is the strongest. As the redshift decreases, heating rate 
fluctuations decrease because of the decreasing bias. 

Figure [4] shows the volume weighted ratio of the in¬ 
ternal energy when blazar heating is included to the in¬ 
ternal energy when blazar heating is not included, as a 
function of the density. For both simulation models, we 
divided the simulation volume into equally-sized small 
sub-volumes and computed the internal energy therein. 
Because both simulations have been started from identi¬ 
cal initial conditions, the interference of the primordial 
density waves gives rise to a comparable morphology of 
the cosmic web so that we can compare thermodynamic 
quantities of (almost) identical sub-volumes of the differ¬ 
ent simulations. These maps clearly highlight that blazar 
heating has more impact in underdense regions, as the 
heating rate per baryon is higher. Even if these regions 
receive less heat than regions with higher density (see 
Figure [3]), blazar heating increases the internal energy by 
a factor two at z = 3 and about an order of magnitude 
at z = 1. At low redshift, even in this inhomogeneous 
model, most of the gas is heated up and there is minimal 
difference between the two models we used for the bias. 

The impact of inhomogeneous heating translates into a 
more complex temperature-density relation, as is shown 
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h- 1 Mpc 





h- 1 Mpc 



Fig. 3.— Slices through the midplane of the box, for z = 3 (left column) and z = 1 (right column): logarithm of the density with respect 
to the mean value (upper row), heating rate fluctuations with the galaxy (middle row) and quasar bias model (lower row). 


on Figure [5j The color map shows the mass weighted 
T — pgas relation from our simulations and the grey con¬ 
tours show the ca se for uniform blazar h eating with the 


same resolution ( 

Puchwein et al. 

2012 

). Temperature 

measures the integrated impact o 

t TeV blazar heating 


over time. The left column shows the T — p gas relation 
in a model with no blazar heating. In the latter case, 
the underdense gas follows a very narrow distribution, 
where the lowest density gas is the coldest. When blazar 
clustering is taken into account, the temperature-density 


relation has a significant scatter for underdense regions 
for z ~ 2 — 3. For the quasar bias model, this scatter 
results in the lower envelope of the temperature-density 
that differs little from the case with no blazar heating. 
However, the mode of the temperature is very close to 
the uniform blazar heating case. At z = 1 both models 
show a very similar behavior and blazar heating can be 
considered nearly homogeneous, though the lower enve¬ 
lope sits at a lower temperature. 

Figure [6] shows the volume-weighted probability dis- 
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Fig. 4.— Ratio between the internal energy (U = 3fcgT/2) in simulations with blazar heating to those without, as a function of overdensity. 
The plots show the galaxy bias model (top) and the quasar bias model (bottom) for z = 3 (left) and z = 1 (right). The color scale is 
logarithmic. 


tribution functions of the temperature for all the simu¬ 
lations. This provides a more quantitative view of the 
scatter in temperature and highlights the impact of clus¬ 
tering on TeV blazar heating. The simulations with in¬ 
homogeneous blazar heating show significant deviation 
from the uniform case for z ^ 2. The higher value of 
the quasar bias results in the presence of more unheated 
gas than in the galaxy bias model. Here the effect of the 
lower envelope is clear as there is a significant tail toward 
lower temperatures: the coldest zones have T < 5 x 10 3 iT 
while the warmest zones have T ~ 10 b K for z = 3. Con¬ 
versely, the warmest gas is only slightly warmer than the 
mode of the temperature. At lower redshift, the mode 
of the temperature is very similar in all models but the 
inhomogeneous models exhibit a larger amount of colder 
gas. 

To have a better understanding of the heating fluctu¬ 
ations with respect to density fluctuations we show the 
mass-weighted Sh — £ ga s distribution on Figure [7| The 
heating rate represents an instantaneous view of The im¬ 
pact of TeV blazar heating as opposed to temperature or 
internal energy which probe the integrated injection his¬ 
tory of non-gravitational heat. As in Figure [3j the quasar 
bias model stands out at high redshift (lower left panel). 
In this model, most of the particles receive slightly more 


heat than in the uniform case. In contrast, the addi¬ 
tional heat is only a few times more than the uniform 
case, while certain regions receive orders of magnitude 
less heat than the mean value. These areas suffer from 
the decrease of the TeV flux and isolation from massive 
structures. This is consistent with the temperature prob¬ 
ability distribution function presenting a prominent low 
temperature tail and a low probability for high temper¬ 
atures for z ^ 2. In the galaxy model, most of the gas 
is heated similarly to the uniform case, translating the 
lower bias for galaxies. 


5. DISCUSSION 

Previous work on TeV blazar heating showed a signifi- 
cant increase in the temperature of the lo w density IGM 
(|Chang et al.|2012||Puchwein et al.|2012|). These results, 
while in good agreement with observational data for the 
mean transmitted flux statistics as well as fits to individ¬ 
ual lines , appear to be in potential tension with recent 
work by Rudie et al. (2012). This observation suggest a 


significant presence of gas in the IGM has not been ex¬ 
posed to additional heating beyond photoheating. In the 
context of blazar heating, this suggests that the heating 
is not entirely uniform. In this work, we have studied the 
impact of inhomogeneous heating and find that when ac- 
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Fig. 5. — Volume-weighted temperature - density relation at z — 4, 3, 2,1 (from top to bottom) for the simulations with no blazar heating 
(left), inhomogeneous heating following galaxy bias (mid dle), and the quasar b ias model (right). The overlying black contours show the 
corresponding T — p gas relation for uniform blazar heating (Puchwein et al.|2012) for the same redshift range. The color scale is logarithmic. 
None of the models includes radiative transfer effects of Hell reiomzation 
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Fig. 6.— Volume-weighted temperature probability distribution 
function for z = 1 (black), 2 (blue), 3 (cyan) and 4 (green) for 
the quasar bias model (solid line), galaxy bias model (dashed line) 
and the uniform case (dotted line). The peaks in the tempera¬ 
ture distribution move towards higher temperatures as the redshift 
decreases. 


counting for the clustering of sources, the variability in 
the heating rate is significant. In particular, we find that 
while most of the gas receives close to average heating, 
and has a temperature close to the uniform heating case, 
there is an important scatter in the temperature of the 
gas for z > 2 . 

The lower envelope of the T — p distribution we find is 
within the error bar s of the powe rlaw fit to the T — p dis¬ 
tribution derived in Bolton et al. (2014). Detailed com¬ 
parison of the Lyman a forest properties will be the sub¬ 
ject of a forthcoming paper. If confirmed with the anal¬ 
ysis of the resulting Lyman a forest, our model would be 
able to account for both the observed increase in temper¬ 


ature in th e low density IGM (Boera et al.||2014 Becker 
et al.|| 20 lT Viel et ah||2009 ) and the lower valu es of the 


line widths as suggested by Rudie et al. (2012). Taken 
together, both observational diagnostics appear to point 
to a more complicated temperature-density relation with 
a possibly broad, multi-valued distribution as opposed to 
a simple power-law relation. 

Around 2.7 < 2 < 4, the impact of TeV blazar heat¬ 
ing is qualitatively similar to the impact of late reion¬ 
ization of Hell. Various physical models and assump¬ 
tions for the quasar distribution can result in highly 
variable temperatu re-density distributions in models of 
Hell reionization (Bolton et al. 2004 |McQuinn et al.f 


2009). However, comparing the temperature-density dis- 
tributions obtained from such models with our Figure [5j 
we notice two major differences. While our model and 
reionization models present a similar lower envelope of 
~ 3 x 10 3 K at p/p = 0.1 (at £ = 3), our model includes a 
significant fraction of gas above a few 10 4 K, which is not 
present in reionization models. On top of this increased 
scatter, our model presents a warmer mean temperature 
than all the Hell models. The same comparison holds 
at z = 4. 

Our model also predicts a warmer IGM at z = 2 , which 
would not be expected from He II reionization, as the de¬ 
pendence on the exact He I I reinozation history will have 
largely faded at that time (Compostella et al. p0l3 ). An 
increased temperature, with a strong scatter m under- 


dense regions could then be attributed to inhomogeneous 
blazar heating. Direct measurements of the low density 
IGM are hard to obtain as the Lyman a forest is sensitive 
to over densities of at least a few. The Lyman a forest 
of He II traces lower density regions and is a promising 


observab le, with more data becoming available (Worseck 


et al. 


Our model neglects important physical aspects of TeV 
blazars such as their duty cycle and beaming of the 
gamma-ray emission. The 7 -ray d uty cycle of B L La c 
objects is of order of 10 % (Stecker & Salamon 1996). 
As TeV blazars have l ow accretio n rates, the ir lifetime is 
long ~ 5 — 7 Gyr ([Cavaliere & D’Elia 2002 ). The 7 ray 
sources are highly b eamed, with observed o pening angles 
peaking around 20° (|Pushkarev et al.|2009 ). Parsec scale 
VLBI observations of the inner parts of radio jets indicate 
variations of the position angle of aro und 20 °, up to 120 ° 
over more than a decade (Lister et al. 2013|). Assuming 
the gamma-ray emission follows the orientation of the 
radio jet, the rapid variability of the sources can poten¬ 
tially increase the scatter in the heating rate. Indeed, in 
the early universe, when few sources are present, indi¬ 
vidual variability increases the inhomogeneity of blazar 
heating. In such case, our model is a lower limit for the 
true scatter. However, later on, as the number of sources 
at a given point increases, the variability of the heating 
rate will be reduced. On top of the intrinsic variability of 
the sources, shot noise, which becomes important when 
sources are rare, is not taken into account and will result 
in additional heating rate fluctuations. For instance, at 
z ~ 1 — 2 , the number of sources that contributes half 
the heating is ~ 10 whereas at z = 3, it dwindles to ~ 1 
(Chang et al. 2012 ). While this would suggest that the 
shot noise is potentially large, this is mitigated somewhat 
by the fact that the next 25% of the heating is provided 
by ~ 100 sources between z = 1 — 2 and ~ 10 source at 
z — 3. Hence, the fluctuations that we calculate in this 
work are likely a lower limit. 

Our model is based on quasar bias, w hich is likely 
lower than TeV-blazar bias (Allevato et al.(2014b. Radio 


lou d AGN are assoc iated with red giant elliptical galax¬ 
ies (|Hyvonen et al. 2007) and are often at the center of 
small clusters and groups. Quasars have a wider distri¬ 
bution of host galaxies and are rarely found at the center 
of clusters. Therefore, the scatter found in our simula¬ 
tions at low redshift is probably a lower limit to the exact 
scatter. In our simulations, any dense source emits TeV 
photons, regardless of the type of galaxy, and history 
of merger, accretion or star formation. However, based 
on the comparison of the galaxy and quasar bias models 
presented here with the uniform heating case we expect 
a more detailed physical model will produce comparable 
results. 


6. CONCLUSIONS 


In this paper we have implemented inhomogeneous 
TeV blazar heating into cosmological simulations to 
study the impact of clustering on the thermal stat e of the 
intergalactfe medium. This extends the work by Chang 


et al. (2012); Puchwein et al. (2012), based on uniform 
blazar heating. 

We developed a filtering function relating heating rate 
fluctuations to the linear dark matter fluctuations. Us¬ 
ing this window function, we are able to model the rele- 





























































11 



Fig. 7.— Volume-weighted distribution of heating rate fluctuations (1 + 5h) with respect to the density fluctuations (1 + S gas ) for the 
galaxy bias model (upper row) and the quasar bias model (lower row) at z = 3 and z = 1. The black dashed line represents the case of 
uniform blazar heating. The color scale is logarithmic. 
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vant length scales for the blazar heating fluctuations and 
include them in a statistical fashion. Our method is a 
cost effective alternative method to fully self consistent 
simulations of blazar heating where modeling black hole 
formation and growth and complete radiative transfer 
would have been prohibitive. 

Our model for the blazar heating fluctuations yields a 
temperature-density relation which can potentially rec¬ 
oncile the observed incre ased mean temperature of the 
IGM (Boera et al. 2014), while maintaining a fraction 
of cold gas, responsib le for the lower envelope of the 
linewidth distribution (Rudie et al. 2012). Therefore, de¬ 
tailed modeling of the Lyman a forest will be the subject 
of a forthcoming paper. If confirmed, this will clearly in¬ 
dicate a more complex thermal history of the IGM, with 
potentially an important impact on late forming struc¬ 
tures. 
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APPENDIX 

We detail the derivation of the window function of Eq. a. To present a transparent derivation of the window 
function (and to clarify some confusion in the literature), we start with a simple toy model to which we add 
progressively more physics. A familiarized reader may directly switch to Appendix [C] In Appendix [Aj we start 
from a purely Newtonian universe and assume that heating rate fluctuations perfectly trace density nuctuations 
and a monochromatic source of very high gamma-rays. Then in Appendix [Bj we include the impact of cosmological 
expansion and a spectral energy distribution of TeV blazars that evolve with redshift as quasars, albeit with a 
significantly smaller normalization in the luminosity density. This impacts on the energy-dependent mean free path 
of TeV photons. Finally, in Appendix [C] we also account for the clustering of sources (which we model through the 
linear bias factor) and redshift distortions as a result of peculiar infall velocities onto cluster potentials. 


A. NEWTONIAN CASE 

A.l. Fluctuations with respect to the mean heating rate 

The TeV flux received (in erg s -1 cm -2 ) at position x, at a given energy, is given by the sum over all the sources 

,2tt roo £( x ') 


pZ7T f*7T f‘OC 

J(x) = / d<J> sinOdO / |x' — x| 2 d|x' — : 

Jo Jo Jo 


4'7r|x / — x|‘ 


= f <m f 

Jn Jo 


£(r + x) _ r/D 

dr - -e r/1Jp 

47r 


(Al) 


where the emissivity £ is an energy per unit time, per unit volume, r = |x' — x|/D pp is the optical depth along the 
line of sight and D pp is the mean free path for pair production, which is constant over the volume. In the last step, 
we introduced r = x' — x, dVt = sin OdOdcf). We have further assumed that emissivity is monochromatic such that: 

&(x) = f(x)5(£7-£7 0 ), (A2) 

where S is the Dirac delta function and Eq is the energy of the monochromatic photons. Similarly the spectral flux is 

Je(x) = J(x)S(E-E 0 ). (A3) 


For simplicity, we have integrated equation (Al) over the monochromatic source. 

The differential heating rate is defined as dQ(x) = dJ(x)/D pp such that the volume-integrated heating rate Q and 
its average (in units of erg cm -3 s -1 ) are given by 

e~ r /Dp 


Q(x) =-!- [ dVt [ dr£( r + x) 
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PP 


The heating rate fluctuations at a given point x are then given by 

<5 H (x) = ^ (X) Z ~ ^ = —4- f dSl f°° dr[£{ r + x) - S]e~ r/D 
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where 5 E are the fluctuations of the TeV photon emissivity. 


(AS) 
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Fig. 8.— Window function for a monochromatic source in a non-expanding universe (Eq. ( |A9| >) . The dashed red line shows a small 
mean free path (D pp = 10“ 1 Mpc), the solid blue line has D pp = 100 Mpc. 


A.2. Window function 

As the universe is infinite and asymptotically flat, we can expand the fluctuations into planar waves, in order to get 
the length scale dependence of heating rate fluctuations, 


i r°° 

1 f‘°° / 

Mr + x) = -—3 / d 3 k'6 E (k')e- ik < r+x \ 
W J- OO 


(A6) 


Performing a convolution of the density fluct uati ons with the kernel Cr 2 e r / Dpp (with C an arbitrary constant), and 
using the Fourier convolution theorem, Eq. (|A5|) rewrites 


~S H (k) = ~S E (k)W(k), (A7) 

where the tilde denotes quantities in Fourier space and W(k) is the Fourier transform of Cr~ 2 e~ r ! Dpp . This yields 


Sh( k) = -—— [ dVt [ dre v ^ Dpp Se( k)e zkr 
47rT)pp Jq J 0 


(AS) 


Introducing fi = cos 6, where is t he angle between the wave vector and the line of sight, and assuming statistical 
isotropy and homogeneity, Eq. (A8) can be simplified as follows, 


5 H (k) = -t- f dn r dre- r ' D ™~5 E {k)e- ik ^ = S E (k)-f- f°° s _^M e ~r/ D 
j-l do ^PP Jo 


? dr 
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y PP J ~ 

S E(k) 


arctan (kD pp ) = S(k) 


arctan (. kD pp ) 
kDr>r> 


In the last step, we have assumed for simplicity that the fluctuations in the TeV emissivity follow the DM density 
fluctuations, i.e., Se = S. While we may justify this by t he tight r elation between the mass of the central black 
hole and the stellar mass of the bulge of the _host galaxy (Haring & Rix 2004), we generalize the mapping between 
emissivity and DM fluctuations in Appendix 0 Figure [8] shows window functions with different mean free paths in a 
purely Newtonian universe. When absorption is present, the impact of large scale structure (i.e. low wave numbers k) 
remains constant as distant sources are absorbed. For scales equal to or larger than the cutoff, heating fluctuations 
follow density fluctuations and overdense regions get more heat. At smaller scales, the density structure has less impact 
on the heating, unless a strong overdensity is present. 

B. EXPANDING UNIVERSE 

The heating rate in an expanding universe at a comoving point x at redshift 2 is given by 


rE m; 
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where the interval E m i n to E ma x is the energy range over which pair production is possible. Quantities indicated by a 
hat are measured in the comoving frame. 

Je is the specific intensity per unit energy (in erg s -1 cm -2 erg -1 ). Assuming a flat geometry on spatial hypersurfaces 
of the Friedmann-Robert son-Walker metric, Je is given by 

Je(E, x) = ^ El J n dn J Z YEz' + ^ z ') e ~ T{E,z,z ’\ (B2) 


47T 


where the vector ?(r(z, z'\ 0, </>), z is the absorption redshift, £e — e#/( 1 + z') 3 is the comoving specific emissivity 
defined in Eq. 0, E' = E(1 + z')/(l + z) and r is the approximate optical depth given by 


t(E,z,z') = t = 


/ dz"- 

J Z J 


dr 


D PP (E, z") dz" 


(B3) 


An alert reader will recognize that the true optical depth must include the effect of a changing E" = E(1 + z")/(l + z) 
as we integrate along z" from z to z'. However, here we assume the typical mean free paths for high energy photons 
are much smaller than the Hubble length so that E" ~ E over the redshift interval, which is small compared to z. We 
then find that the energy-integrated and volume-integrated heating rate is then given by 

C[*M1 = r /«/” (B4) 
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where 7t = x + r. Similarly to Eq. ( |A4| ), the mean heating rate is given by 
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This enables us to define the heating rate fluctuations 




Q[St(z)] - Q(z) (1 + z) 3 c 


f 

Je x 

p r~ „ Z \ [ / 

3 JE mir , Dpp{E,z) J n J z 


dE 


Q(z) 

(1 +z) 3 c /,Bmax 

47 tQ 


47 tQ 
dE 


D V p(E , z) 


[ dO, [ 

Jn Jz 


dz' 


S E >(1Z,z')-£ e ,(z') 


e -T( E ,Z, Z ') 


(1 + z')H{z') 


r0 ° dz , £E'{z')5 E {ii,z')e-^ E ^') 


(1 + zf) H{zJ ) 


(B6) 


The TeV emission is related to the presence of supermassive black holes at the center of galaxies, which are located 
in collapsed dark matter halos. We can thus connect the fluctuations of the TeV emission, within a certain radius f, 
to the underlying dark matter fluctuations 5. At this stage we assume that TeV fluctuations exactly match the dark 
matter fluctuations 5 E = 6. We explain in the next section that this is not exactly true and will take into account 
various corrections. The initial density fluctuations represent a Gaussian r andom field, wh ose exa ct properties depend 
on the earliest stages of the universe pri or to recombination ( Bardeen et al.|[l986 Peebles 1982). They grow linearly 
between z' and z following (Heath|[l977) 


5('Tt,z)=S(n,z)D(z)/D(z'), 


where the linear growth factor is given by 


d W = d 0 *( 2) /“*'1±T 


(B7) 


(B8) 


The linear approximation breaks down when the amplitude of the root mean square of th e pert urbati ons ap proaches 
unity. The evolution of the density field is then determined by the spherical collapse (Gunn & Gott 1972) and the 
virialization of halos. In the linear regime, the growth of the modes is independent of the wavenumber and we have 

D(z') D(z') 


S E (r, z') = (5(r, z') = 5( r, z) 


D{z) D{z) (2 


ip / <f>k'4(k')< 


/k'r 


And Eq. (B6) rewrites to 
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(B9) 


(BIO) 


Four ier t ransforming yields <5(k) on the left hand side while the right hand side transforms in a similar fashion to 
Eq. (A9). As the power spectrum of density fluctuations is statistically isotropic and homogeneous, this simplifies to 


x n \ 1ft + z ) 3c T max dE [ 

Sa(k , z) = s(k , e) —^j^ 
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Figure [9] shows the window function in an expanding universe at different redshifts. 


The position of the breaks 
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Fig. 9.— Window function for blazar heating in an expanding universe without account for clustering bias (Eq. ( |B11| >) for z — 1 (solid 
line), 2 (dashed line) and 4 (dotted line). 


results from a combination between the value of the mean free path and the redshift evolution of the blazar luminosity 
density. As discussed in the main text, the rapidly decreasing blazar luminosity density with increasing redshift for 
z > 2, which moves the break toward larger wave vectors, dominates over the increasing comoving mean free path, 
which moves the break toward smaller wave vectors. This is due to the fact that the integration over the energy of the 
photons, which includes lower energy photons with larger mean free paths, allows the former effect to dominate and 
moves the break toward smaller wave vectors with increasing redshift. The TeV emission fluctuations are not exactly 
equal to the dark matter density fluctuations. In the next section we will account for the various corrections that have 
to be taken into account to determine a more accurate window function. 


C. COMPLETE WINDOW FUNCTION 

The TeV fluctuations are biased with respect to the dark matter fluctuations, as we detail in §2.2[ yielding 

S E (Z,z) = b(z)6(k,z), 


(Cl) 


with b the Eulerian bias. If galaxies were moving exactly with the Hubble flow, their redshift would yield their exact 
distance to an observer. However, structures falling towards a central potential of a cluster have an infall velocity 

towards the central overdensity, yielding different overdensities in redshift space 6$ and real space. In our case, the 
blazar luminosity density evolution, the mean free path, and the bias as a function of redshift are all subject to this 
redshift space distortions as measu red by_the z = 0 observer. We follow the derivation of the redshift space distortions 
of density fluctuations of Mo et al. (2011). The key difference that we encounter is that redshift space distortions 
impact on H = x + ?, whereas our computation involves the displacements ? between the heated point and a distant 
source. To this end, we define a comoving redshift space distance 


A - ^ 

S = aH(z) 


= a + 


Viz 

H(z) 


&K, 


(C2) 


which is given in units of length, a = 1/(1 + z) is the cosmological scale factor, v E is the radial component of the 
comoving peculiar velocity (v = v/a with v denoting the physical peculiar velocity), and e E is the unit vector along 
the line of sight. The conservation of TeV blazar number then implies 


hg\s)d 3 s — n E {H)d 3 H, or 


i+4 s) (§) 


i 3 8 = 


1 + 


few] 


d 3 n , 


(C3) 


where we employed the fact that the mean number density is the same in both spaces. Assuming statistical isotropy, 
we find 


i+4°(s) = 


i + 


few] 
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( H + vn/H(z)y 
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dvn \ 
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Assuming \vn\ <C 7 ZH(z) and keeping only terms up to lowest non-trivial order, we obtain 


[§(£)] = s E (n) - 


d 

^ ) vn ~ ^e(TZ) — 


h{z) \dn n 


1 dv n 


(C5) 


where we assumed in the second step that the scale of the fluctuations ks small compared to the distance to the sources 
(i.e., we used the plane-parallel approximation according to [Kaiser 1987). In the li near regi me (S E ( IV) <C 1), the 


peculiar velocity field is related to the density field via the Zerdovich approximation ( |Mo et a l 

- f -v b dw 

5e(71) = -- 


me [o E [ 

l2onJ, 


fH(z)dH’ (C6) 

where / = d log S/d log a. It can be shown in the linear regime that non- trivi al solutions of the velocity field must be 
irrotational so that v = Using this property, we can rewrite Eq. ( |C 6 | ) to yield 

fH{z) d 


vn = -- 




where V ~ 2 represents the inverse Laplacian. It follows that 


iz 


4* ) [^(f)] = 


1 


f-?Lv^ 2 1 S E [H( f)] = 
bdK 2 
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lAv - 2 

6 <9f 2 f 


5 e [K(v)}. 


(C7) 


(C 8 ) 


This demonstrates that redshift space distortions do not depend on the line-of-sight distance from the observer to 
the source but only on the shape of the local gravitational potential of the source! This property enables us to treat 
redshift space distortions locally by performing a Fourier decomposition of the density field into plane waves, which 
yields 


Ks) 


k (J) (k) = ( 1 + 0 2 ) fe(k) = (b + j> 2 ) <5(k), 


where /i = k z /k. Keeping only first order corrections for density fluctuations, Eq. (p4j) yields 




dz' 


47r JE min D pp (E,z)J n J z (l+z')H(z') 

Subtracting the mean heating rate (Eq. (|B5|)) then yields the fluctuations 




where we introduced 
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(C 12 ) 


for convenience. Transforming to Fourier space yields 

,X(E, z, z') 


1 ^max p p 

S H {k,z) = —=■ / dE dSl 

4:7V (D J ^min 


dz 
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Assuming statistical isotropy and homogeneity allows us to simplify this expression, and we arrive at 

S H (z,k) = ^k_ f max dE f dfl [ dz'X{E,z,z')^Ed f dfi \ e ~ ik ^(b + f/j, 2 ) 
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where we used the spherical Bessel functions 
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The window function for the heating rate fluctuations is then given by 
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(C19) 


This equation is similar to |Pritchard & Furlanetto| (|2007|) and |Barkana & Loeb| (2005) 
unnecessary area correction factor, which would be present only for a Lagrangian aescripti 


but does not contain the 
Lagrangian description of bias. 
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